
    DeltaEst = function(x.delta, x.delta.d, x.op.y, x.op.d, x.pi, z, d, y, 
                        weights = NULL){
        # x.delta = x.delta.d = x.op.y = x.op.d = x.pi = X; z=Z; d=D; y=Y;  weights = wt
        
        if(is.null(weights)) weights = rep(1,length(y))
        
        ptm = proc.time()
        
        Delta.est = numeric(length(method))
        names(Delta.est) = method
        
        ### 1 = regression estimator
        system.time({reg.sol = Reg(x.delta, x.delta.d, x.op.y, x.op.d, z, d, y,
                                   weights)})
        Delta.est["reg"] = reg.sol$Delta.est
        
        ### 2 = IPW estimator
        system.time({ipw.sol = IPW(x.pi, x.delta.d, x.delta, z, d, y, weights)})
        Delta.est["ipw"]   = ipw.sol$Delta.est
        Delta.est["b-ipw"] = ipw.sol$b.Delta.est
        
        ### 3 = G-estimator
        system.time({g.sol = GEst(x.delta, x.pi, z,d,y, weights)})
        Delta.est["g"] = g.sol$Delta.est
        
        ### 4 = Triply robust estimator
        system.time({tr.sol = TR(x.delta, x.delta.d, x.op.y, x.op.d, x.pi, z, d, 
                                 y, reg.sol, ipw.sol, weights)})
        Delta.est["tr"]   = tr.sol$Delta.est
        Delta.est["b-tr"] = tr.sol$b.Delta.est
        
        time = (proc.time() - ptm)[3]
        
        return(list(Delta.est=Delta.est, reg.sol=reg.sol, ipw.sol=ipw.sol,
                    g.sol=g.sol, tr.sol=tr.sol, time = time))
        
    }